library(data.table)
library(plyr)
library(dplyr)
library(ggplot2)
library(foreign)
library(doBy)
library(magrittr)
library(gmodels)
library(survey)
library(tidyr)
library(RColorBrewer)
library(ggthemes)

####################################################################################################################
##################List of Included Files############################################################################
####################################################################################################################
# eckhouse-2019-jop-analysis.R: data analysis file

#GSS data
#gss <- read.csv(file = "gss4jop.csv")
#codebook available on the GSS website: gss.norc.org


#AAMS data
#afam <- read.csv(file = "aams4jop.csv")
#Codebook: AAMS Codebook.pdf


# AP-NORC data 
# norc <- read.csv("norc4jop.csv", stringsAsFactors = F, na.strings = c("-99", ""))
# Codebook: "APNORC Public Use Files Codebook.pdf"




####################################################################################################################
##################GSS DATA ANALYSIS (FIGURES 1 & 2, APPENDIX A)#####################################################
####################################################################################################################
gss <- read.csv(file = "gss4jop.csv")


#data processing 
gss$RACEcode <- mapvalues(gss$RACE, from = c(1, 2, 3), to = c("white", "black", "other"))
gss$RACEcode <- as.factor(gss$RACEcode)
gss$POLHITOKcode <- as.numeric(mapvalues(gss$POLHITOK, from = c(1, 2, 8, 9, 0), to = c(1, 0, NA, NA, NA)))
gss$PARTYIDlean <- mapvalues(gss$PARTYID, from = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), 
                             to = c("Democrat", "Democrat", "Lean Dem", "Independent", "Lean Rep",  "Republican", "Republican", "Other", "Don't know", NA))	

gss$PARTYIDlean <- mapvalues(gss$PARTYID, from = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), 
                             to = c("Democrat", "Democrat", "Democrat", "Independent", "Republican",  "Republican", "Republican", "Other", "Don't know", NA))	

gss$PARTYIDlean <- as.factor(gss$PARTYIDlean)
#courts question: recoded to, basically 'should the courts be more harsh (1), less harsh (-1), or about right (0)
gss$COURTScode <- as.numeric(mapvalues(gss$COURTS, from = c(1, 2, 3, 8, 9, 0), to = c(-1, 1, 0, NA, NA, NA)))

#capital punishment
gss$CAPPUNcode <-as.numeric(mapvalues(gss$CAPPUN, from = c(1,2,8,9,0), to = c(1, 0 , NA, NA, NA)))

sort(grep("NAT", names(gss), value = T))


#### SPENDING QUESTIONS ######
#foreign aid, should we spend more (1), less (-1), same (0)
gss$NATAIDYcode <- mapvalues(gss$NATAIDY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATAIDcode <- mapvalues(gss$NATAIDY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#Military/armaments/defense
gss$NATARMSYcode <- mapvalues(gss$NATARMSY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATARMScode <- mapvalues(gss$NATARMS, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#assistance for childcare
gss$NATCHLDcode <- mapvalues(gss$NATCHLD, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#solving the problems of big cities
gss$NATCITYcode <- mapvalues(gss$NATCITY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATCITYYcode <- mapvalues(gss$NATCITYY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
#halting the crime rate: should we spend more
gss$NATCRIMEcode <- mapvalues(gss$NATCRIME, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
#law enforcement
gss$NATCRIMYcode <- mapvalues(gss$NATCRIMY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
#dealing with drug addiction
gss$NATDRUGcode <- mapvalues(gss$NATDRUG, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATDRUGYcode <- mapvalues(gss$NATDRUGY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

gss$NATEDUCcode <- mapvalues(gss$NATEDUC, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATEDUCYcode <- mapvalues(gss$NATEDUCY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATENRGYcode <- mapvalues(gss$NATENRGY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

gss$NATENVIRcode <- mapvalues(gss$NATENVIR, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATENVIYcode <- mapvalues(gss$NATENVIY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#welfare
gss$NATFAREcode <- mapvalues(gss$NATFARE, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATFAREYcode <- mapvalues(gss$NATFAREY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#health care, pretty hilarious change post-Obamacare
gss$NATHEALcode <- mapvalues(gss$NATHEAL, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATHEALYcode <- mapvalues(gss$NATHEALY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#mass transit
gss$NATMASScode <- mapvalues(gss$NATMASS, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#parks
gss$NATPARKcode <- mapvalues(gss$NATPARK, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
#improving conditions of blacks
gss$NATRACEcode <- mapvalues(gss$NATRACE, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATRACEYcode <- mapvalues(gss$NATRACEY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#roads
gss$NATROADcode <- mapvalues(gss$NATROAD, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#social security
gss$NATSOCcode <- mapvalues(gss$NATSOC, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#science
gss$NATSCIcode <- mapvalues(gss$NATSOC, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))

#space
gss$NATSPACcode <- mapvalues(gss$NATSPAC, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))
gss$NATSPACYcode <- mapvalues(gss$NATSPACY, from = c(1, 2, 3, 8, 9, 0), to = c(1, 0, -1, NA, NA, NA))


#support for affirmative action
gss$AFFRMACTcode <- mapvalues(gss$AFFRMACT, from = c(1, 2, 3, 4, 8, 9, 0), to = c(4, 3, 2, 1, NA, NA, NA))

gss$CAPPUNcode <- mapvalues(gss$CAPPUN, from = c(1, 2, 8, 9, 0), to = c(1, 0, NA, NA, NA))

#MAKE summaries

qsum <- gss %>% filter(!is.na(gss$raceparty)) %>% group_by(YEAR, raceparty) %>% 
  summarise_at(vars(qcoded.y), funs(weighted.mean(., w = WTSSALL), sd, std.error), na.rm=T)



qsum.long <- melt(qsum, id = c("YEAR", "raceparty"))


natmeans <- unique(grep("^NAT.+mean", qsum.long$variable, value = T))

nats <- sort(grep("NAT.+code", names(gss), value = T))


neutralvars <- c("NATAIDcode", "NATAIDYcode", "NATARMScode", "NATARMSYcode", "NATCHLDcode", "NATEDUCcode",  "NATEDUCYcode",
                 "NATENRGYcode", "NATENVIRcode", "NATENVIYcode","NATHEALcode","NATHEALYcode", "NATMASScode", 
                 "NATPARKcode", "NATROADcode", "NATSCIcode",
                 "NATSOCcode", "NATSPACcode", "NATSPACYcode")


crimevars <- c("NATCRIMEcode", "NATCRIMYcode")
racevars <- c("NATRACEcode", "NATRACEYcode")

natslisted <- c(neutralvars, racializedvars, crimevars, racevars)
setdiff(nats, natslisted)

racemeans <- gsub("$", "_weighted.mean", racevars)
racelabels <- c("NATRACEcode_weighted.mean" = "Improve Conditions \n of Blacks",
                "NATRACEYcode_weighted.mean" = "Assistance to Blacks")


neutralmeans <- gsub("$", "_weighted.mean", neutralvars)
neutralmeans <- neutralmeans[!neutralmeans %in% grep("Y", neutralmeans, value = T)]

neutrallabels <- c("NATAIDcode_weighted.mean" = "Foreign aid",  "NATARMScode_weighted.mean" = "Military", 
                   "NATCHLDcode_weighted.mean" = "Childcare assistance",
                   "NATEDUCcode_weighted.mean" = "Education", "NATENVIRcode_weighted.mean" = "Environment", 
                   "NATHEALcode_weighted.mean" = "Health care", "NATMASScode_weighted.mean" = "Mass transit",
                   "NATPARKcode_weighted.mean"  = "Parks & Recreation", "NATROADcode_weighted.mean" = "Roads",
                   "NATSCIcode_weighted.mean" = "Science", "NATSOCcode_weighted.mean" = "Social Security",
                   "NATSPACcode_weighted.mean" = "Space")


racializedvars <- c("NATCITYcode", "NATCITYYcode", "NATDRUGcode", "NATDRUGYcode", "NATFAREcode")
racializedmeans <- gsub("$", "_weighted.mean", racializedvars)
racializedlabels <- c("NATCITYcode_weighted.mean" = "Solving problems \n of big cities", "NATCITYYcode_weighted.mean" = "Assistance to \n big cities", 
                      "NATDRUGcode_weighted.mean" = "Dealing with \n drug addiction", "NATDRUGYcode_weighted.mean" = "Drug rehabilitation",
                      "NATFAREcode_weighted.mean" = "Welfare" )

crimemeans <- gsub("$", "_weighted.mean", crimevars)
crimelabels <- c("NATCRIMEcode_weighted.mean" = "Halting rising \n crime rate", "NATCRIMYcode_weighted.mean" = "Law Enforcement")

otherracemeans <- c("AFFRMACTcode_weighted.mean")
otherracelabels <- c("AFFRMACTcode_weighted.mean" = "Support for \n affirmative action")

othercjmeans <- c("CAPPUNcode_weighted.mean", "POLHITOKcode_weighted.mean", "COURTScode_weighted.mean")
othercjlabels <- c("CAPPUNcode_weighted.mean" = "Support for \n capital punishment", "POLHITOKcode_weighted.mean" = "Concern about \n police violence", "COURTScode_weighted.mean" = "Courts should be \n more/less harsh")


welfarevars <- c("NATFAREcode", "NATCHLDcode", "NATHEALcode", "NATSOCcode")
welfaremeans <- gsub("$", "_weighted.mean", welfarevars)

varlabels <- c(crimelabels, otherracelabels, othercjlabels, racializedlabels, racelabels, neutrallabels)


####Figure 1
ggplot(data = qsum.long[qsum.long$variable %in% c(crimemeans, "YEAR", "raceparty"), ], aes(x = YEAR, y = value)) + 
  geom_point(aes(shape = raceparty), size = 3) + geom_line(aes(linetype = raceparty)) + 
  facet_wrap(~variable, labeller=as_labeller(varlabels)) + theme_tufte(20) +
  coord_cartesian(ylim = c(0.25,1), xlim = c(1974, 2016)) + xlab("Year between 1974 and 2016") + ylab("More or less money to....") + 
  theme_tufte(base_size =15) + scale_shape_discrete(name = "Race & Party") + scale_linetype_discrete(name = "Race & Party")
ggsave("natcrime_x_year_by_raceparty.pdf", scale = .8, width = 12)


####Figure 2
ggplot(data = qsum.long[qsum.long$variable %in% c(othercjmeans, "YEAR", "raceparty"), ], aes(x = YEAR, y = value)) + 
  geom_point(aes(shape = raceparty), size = 2) + geom_line(aes(linetype = raceparty)) + 
  facet_wrap(~variable, ncol = 3, scales = "fixed", labeller = as_labeller(othercjlabels)) + 
  coord_cartesian(xlim = c(1974, 2016)) + xlab("Year") + ylab("Criminal Justice Questions") + 
  theme_tufte(base_size =15) + scale_shape_discrete(name = "Race & Party") + scale_linetype_discrete(name = "Race & Party")
ggsave("othercj_x_year_by_raceparty.pdf", scale = .8, width = 18)

####Figure 5 (Appendix A)
ggplot(data = qsum.long[qsum.long$variable %in% c(neutralmeans, "YEAR", "raceparty"), ], aes(x = YEAR, y = value)) + 
  geom_point(aes(shape = raceparty), size = 2) + geom_line(aes(linetype = raceparty)) + 
  facet_wrap(~variable, ncol = 3, labeller = as_labeller(c(neutrallabels, racializedlabels)), scales = "free_y") + theme_tufte(20) + 
  coord_cartesian(xlim = c(1974, 2016)) + xlab("Year between 1974 and 2016") + ylab("More or less money to....") + 
  theme_tufte(base_size =15) + scale_shape_discrete(name = "Race & Party") + scale_linetype_discrete(name = "Race & Party")
ggsave("natneutral_x_year_by_raceparty.pdf", scale = .8, width = 18, height = 18)

####Figure 6 (Appendix A)
ggplot(data = qsum.long[qsum.long$variable %in% c(racializedmeans, racemeans, "YEAR", "raceparty"), ], aes(x = YEAR, y = value)) + 
  geom_point(aes(shape = raceparty), size = 2) + geom_line(aes(linetype = raceparty)) + 
  facet_wrap(~variable, ncol = 3, labeller = as_labeller(c(racializedlabels, racelabels)), scales = "free_y") + theme_tufte(20) + 
  coord_cartesian(xlim = c(1974, 2016)) + xlab("Year between 1974 and 2016") + ylab("More or less money to....") + 
  theme_tufte(base_size =15) + scale_shape_discrete(name = "Race & Party") + scale_linetype_discrete(name = "Race & Party")
ggsave("natracializedrace_x_year_by_raceparty.pdf", scale = .8, width = 18, height = 18)






####################################################################################################################
##################AFRICAN AMERICAN MEN'S SURVEY DATA ANALYSIS (FIGURE 3, TABLE 1)###################################
####################################################################################################################
afam <- read.csv("aams4jop.csv")

#data processing
#QS4, race
racebegin <- levels(afam$QS4)
afam$race <- afam$QS4

afam$race <- mapvalues(afam$race, from = racebegin, to = c("White", "Black", "Asian", "Other", "Don't know", "Refused"))
rm(racebegin)

#53 why more black men in jail, is X a big reason
reason <- levels(afam$Q53A)
reasonint <- c(2, 1, 0, NA, NA)
afam$Q53Aint <- as.numeric(as.character((mapvalues(afam$Q53A, from = reason, to = reasonint))))
afam$Q53Bint <- as.numeric(as.character((mapvalues(afam$Q53B, from = reason, to = reasonint))))
afam$Q53Cint <- as.numeric(as.character((mapvalues(afam$Q53C, from = reason, to = reasonint))))
afam$Q53Dint <- as.numeric(as.character((mapvalues(afam$Q53D, from = reason, to = reasonint))))
afam$Q53Eint <- as.numeric(as.character((mapvalues(afam$Q53E, from = reason, to = reasonint))))
afam$Q53Fint <- as.numeric(as.character((mapvalues(afam$Q53F, from = reason, to = reasonint))))
afam$Q53Gint <- as.numeric(as.character((mapvalues(afam$Q53G, from = reason, to = reasonint))))

#personal responsibility: b & e
afam$incarcerationpersonal <- rowSums(subset(afam, select = c(Q53Bint, Q53Eint)), na.rm = T)/4
afam$incarcerationsocial <- rowSums(subset(afam, select = c(Q53Aint, Q53Fint, Q53Gint)), na.rm = T)/6
afam$incarcerationbias <- rowSums(subset(afam, select = c(Q53Cint, Q53Dint)), na.rm = T)/4

#Party ID
partystrong <- levels(afam$Q70)
afam$partyIDstrong <- mapvalues(afam$Q70, from = partystrong, to = c("Republican", "Democrat", "Independent", "Other", NA, NA))
afam$partyIDlean <- mapvalues(afam$Q70, from = partystrong, to = c("Republican", "Democrat", NA, NA, NA, NA))

aasvy <- svydesign(ids = ~0, data = afam, weights = afam$WEIGHT)

#TABLE 1

svymean(incarcerationbias~0, subset(aasvy, race == "Black" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(incarcerationbias~0, subset(aasvy, race == "White" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(incarcerationbias~0, subset(aasvy, race == "White" & partyIDlean == "Republican"), na.rm = TRUE)


svymean(incarcerationsocial~0, subset(aasvy, race == "Black" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(incarcerationsocial~0, subset(aasvy, race == "White" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(incarcerationsocial~0, subset(aasvy, race == "White" & partyIDlean == "Republican"), na.rm = TRUE)


svymean(incarcerationpersonal~0, subset(aasvy, race == "Black" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(incarcerationpersonal~0, subset(aasvy, race == "White" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(incarcerationpersonal~0, subset(aasvy, race == "White" & partyIDlean == "Republican"), na.rm = TRUE)



#GRAPHS FOR FIGURE 3

xpersonalybias <- svyby(~incarcerationbias, by = ~incarcerationpersonal + raceparty, design = subset(aasvy, raceparty %in% racepartylabels), FUN = svymean, na.rm = TRUE)
ggplot(data = xpersonalybias, 
       aes(x = incarcerationpersonal, y = incarcerationbias, 
           group = interaction(incarcerationpersonal, raceparty))) +
  geom_point(aes(shape = raceparty), size = 8) +
  geom_errorbar(aes(ymin=incarcerationbias - 1.96*se, ymax=incarcerationbias + 1.96*se), width = 0) +
  xlim(0, 1) + ylim(0, 1) + theme_tufte(20) + 
  scale_shape_manual(name = "Race & Party", labels = racepartylabels,
                     breaks = race.party, values = race.partyshapes) + 
  xlab("Personal Responsibility \n & Incarceration Index") +
  ylab("Racial Bias \n & Incarceration Index")
ggsave("../output/xpersonalybias.pdf")

getwd()
#personal vs social 
xpersonalysocial <- svyby(~incarcerationsocial, by = ~incarcerationpersonal + raceparty, design = subset(aasvy, raceparty %in% racepartylabels), FUN = svymean, na.rm = TRUE)

ggplot(data = xpersonalysocial, 
       aes(x = incarcerationpersonal, y = incarcerationsocial, 
           group = interaction(incarcerationpersonal, raceparty))) +
  geom_point(aes(shape = raceparty), size = 8) +
  geom_errorbar(aes(ymin=incarcerationsocial - 1.96*se, ymax=incarcerationsocial + 1.96*se), width = 0) +
  xlim(0, 1) + ylim(0, 1) + theme_tufte(20) + 
  scale_shape_manual(name = "Race & Party", labels = racepartylabels,
                     breaks = race.party, values = race.partyshapes) + 
  xlab("Personal Responsibility \n & Incarceration Index") +
  ylab("Socioeconomic Factors \n & Incarceration Index")
ggsave("../output/xpersonalysocial.pdf")


xsocialybias <- svyby(~incarcerationbias, by = ~incarcerationsocial + raceparty, design = subset(aasvy, raceparty %in% racepartylabels), FUN = svymean, na.rm = TRUE)
ggplot(data = xsocialybias, 
       aes(x = incarcerationsocial, y = incarcerationbias, 
           group = interaction(incarcerationsocial, raceparty))) +
  geom_point(aes(shape = raceparty), size = 8) +
  geom_errorbar(aes(ymin=incarcerationbias - 1.96*se, ymax=incarcerationbias + 1.96*se), width = 0) +
  xlim(0, 1) + ylim(0, 1) + theme_tufte(20) + 
  scale_shape_manual(name = "Race & Party", labels = racepartylabels,
                     breaks = race.party, values = race.partyshapes) + 
  xlab("Socioeconomic Factors \n & Incarceration Index") +
  ylab("Racial Bias \n & Incarceration Index")
ggsave("../output/xsocialybias.pdf")


####################################################################################################################
##################AP-NORC DATA ANALYSIS (FIGURE 4)##################################################################
####################################################################################################################
#load data set from csv
norc <- read.csv("norc4jop.csv", stringsAsFactors = F, na.strings = c("-99", ""))
head(norc)


#recode relevant variables


#Q33 and Q34 ask about specific policy prefs  
effective <- c(sort(unique(norc$Q33A)), NA)
effint <- c(5, 4, 3, 2, 1, NA)

norc$Q33Aint <-  mapvalues(norc$Q33A, from = effective, to = effint)
norc$Q33Bint <-  mapvalues(norc$Q33B, from = effective, to = effint)
norc$Q33Cint <-  mapvalues(norc$Q33C, from = effective, to = effint)
norc$Q33Dint <-  mapvalues(norc$Q33D, from = effective, to = effint)


norc$Q34Aint <-  mapvalues(norc$Q34A, from = effective, to = effint)
norc$Q34Bint <-  mapvalues(norc$Q34B, from = effective, to = effint)
norc$Q34Cint <-  mapvalues(norc$Q34C, from = effective, to = effint)
norc$Q34Dint <-  mapvalues(norc$Q34D, from = effective, to = effint)
norc$Q34Eint <-  mapvalues(norc$Q34E, from = effective, to = effint)

#Q35 reasons for police violence
reason <- unique(norc$Q35A)
reasint <- c(0, 1, NA, 2)

norc$Q35Aint <-  mapvalues(norc$Q35A, from = reason, to = reasint)
norc$Q35Bint <-  mapvalues(norc$Q35B, from = reason, to = reasint)
norc$Q35Cint <-  mapvalues(norc$Q35C, from = reason, to = reasint)
norc$Q35Dint <-  mapvalues(norc$Q35D, from = reason, to = reasint)
norc$Q35Eint <-  mapvalues(norc$Q35E, from = reason, to = reasint)
norc$Q35Fint <-  mapvalues(norc$Q35F, from = reason, to = reasint)

#make the integer ones numeric

norc <- norc %>% mutate_each(funs(as.numeric), ends_with("int"))



#party + indep

#raceth
racefirst <- unique(norc$raceth)
norc$race <- mapvalues(norc$raceth, from = racefirst, to = c("White", "Black", "Other", "Hispanic"))

party1 <- unique(norc$party)



#mapvalues to D/R/I
norc$partyIDlean <- mapvalues(norc$party, from = party1, to = c("I", "I", "Democrat", "Republican", NA))

#use indx to select I entries only.
indx <- norc$partyIDlean == "I" & !is.na(norc$partyIDlean)
summary(indx)

lean <- unique(norc$indep)
#recode by lean

norc$partyIDlean[indx] <- mapvalues(norc$indep[indx], from = lean, to = c("Republican", NA, NA, "Democrat"))
head(cbind(norc$party, norc$indep, norc$partyIDlean))

norc$raceparty <- paste(norc$race, norc$partyIDlean, sep = " ")
unique(norc$raceparty)
indx <- !(norc$raceparty %in% c("White Republican", "White Democrat", "Black Democrat"))
norc$raceparty[indx] <- NA

#make survey design object
norcsvy <- svydesign(ids = ~0, data = norc, weights = norc$WEIGHT2)
unique(norc$raceth)

svymean(Q33A~0, subset(norcsvy, race == "White", partyIDlean == "Republican"), na.rm = TRUE)
svymean(Q33A~0, subset(norcsvy, race == "White" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(Q33A~0, subset(norcsvy, race == "Black" & partyIDlean == "Democrat"), na.rm = TRUE)


svymean(Q33B~0, subset(norcsvy, race == "White", partyIDlean == "Republican"), na.rm = TRUE)
svymean(Q33B~0, subset(norcsvy, race == "White" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(Q33B~0, subset(norcsvy, race == "Black" & partyIDlean == "Democrat"), na.rm = TRUE)


svymean(Q33Bint~0, subset(norcsvy, race == "White", partyIDlean == "Republican"), na.rm = TRUE)
svymean(Q33Bint~0, subset(norcsvy, race == "White" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(Q33Bint~0, subset(norcsvy, race == "Black" & partyIDlean == "Democrat"), na.rm = TRUE)


svymean(Q33C~0, subset(norcsvy, race == "White", partyIDlean == "Republican"), na.rm = TRUE)
svymean(Q33C~0, subset(norcsvy, race == "White" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(Q33C~0, subset(norcsvy, race == "Black" & partyIDlean == "Democrat"), na.rm = TRUE)


svymean(Q33Cint~0, subset(norcsvy, race == "White", partyIDlean == "Republican"), na.rm = TRUE)
svymean(Q33Cint~0, subset(norcsvy, race == "White" & partyIDlean == "Democrat"), na.rm = TRUE)
svymean(Q33Cint~0, subset(norcsvy, race == "Black" & partyIDlean == "Democrat"), na.rm = TRUE)


svymean(Q33B~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q33B~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q33B~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)


svymean(Q33C~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q33C~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q33C~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)

svymean(Q33D~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q33D~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q33D~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)


svymean(Q34A~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q34A~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q34A~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)

svymean(Q34B~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q34B~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q34B~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)

cor(norc$Q34Aint[norc$race == "Black"], norc$Q34Bint[norc$race == "Black"], use = "complete.obs")
cor(norc$Q34Aint[norc$race == "White" & norc$partyIDlean == "Republican"], 
    norc$Q34Bint[norc$race == "White" & norc$partyIDlean == "Republican"], use = "complete.obs")
cor(norc$Q34Aint[norc$race == "White" & norc$partyIDlean == "Democrat"], 
    norc$Q34Bint[norc$race == "White" & norc$partyIDlean == "Democrat"], use = "complete.obs")

cor(norc$Q34Bint[norc$race == "Black" & norc$partyIDlean == "Democrat"], 
    norc$Q34Cint[norc$race == "Black" & norc$partyIDlean == "Democrat"], use = "pairwise.complete.obs")
cor(norc$Q34Bint[norc$race == "White" & norc$partyIDlean == "Democrat"], 
    norc$Q34Cint[norc$race == "White" & norc$partyIDlean == "Democrat"], use = "complete.obs")
cor(norc$Q34Bint[norc$race == "White" & norc$partyIDlean == "Republican"], 
    norc$Q34Cint[norc$race == "White" & norc$partyIDlean == "Republican"], use = "complete.obs")

cor(norc$Q34Cint[norc$race == "Black" & norc$partyIDlean == "Democrat"], 
    norc$Q34Dint[norc$race == "Black" & norc$partyIDlean == "Democrat"], use = "complete.obs")
cor(norc$Q34Cint[norc$race == "White" & norc$partyIDlean == "Democrat"], 
    norc$Q34Dint[norc$race == "White" & norc$partyIDlean == "Democrat"], use = "complete.obs")
cor(norc$Q34Cint[norc$race == "White" & norc$partyIDlean == "Republican"], 
    norc$Q34Dint[norc$race == "White" & norc$partyIDlean == "Republican"], use = "complete.obs")



svymean(Q34C~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q34C~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q34C~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)

svymean(Q35A~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q35A~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q35A~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)

svymean(Q35B~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q35B~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q35B~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)

svymean(Q35C~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q35C~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q35C~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)

svymean(Q35D~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q35D~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q35D~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)

svymean(Q35E~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q35E~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q35E~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)

svymean(Q35F~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(2) Republican"), na.rm = TRUE)
svymean(Q35F~0, subset(norcsvy, raceth == "(1) White, non-Hispanic" & party == "(1) Democrat"), na.rm = TRUE)
svymean(Q35F~0, subset(norcsvy, raceth == "(2) Black, non-Hispanic"), na.rm = TRUE)




ggplot(data = subset(norc, race %in% c("White") & !is.na(partyIDlean)), aes(x = Q48Bint, y = Q48Dint, group = interaction(Q48Bint, race, partyIDlean))) +
  stat_summary(fun.data = mean_cl_normal, aes(color = interaction(race, partyIDlean))) +
  stat_summary(fun.data = mean_cl_normal, data = subset(norc, race %in% c("Black")), aes(color = race, group = Q48Bint)) +  
  xlim(0, 2) + ylim(0, 2) + theme_tufte()


#######Q33 and Q34 are specific policies
#make tables
eff <- sort(unique(norc$Q34A))
efflevels <- c("Not at all effective", "Not too effective", "Somewhat effective", 
               "Very effective", "Extremely effective")

pal <- brewer.pal(5, "RdBu")
pal <- c(pal[1], pal[5], pal[4], "#DFDFDF", pal[2])

pal <- brewer.pal(5, "Greys")
pal <- c(pal[1], pal[5], pal[4], pal[3], pal[2])

lows <- c("Somewhat effective", "Not too effective", "Not at all effective")
highs <- c("Somewhat effective", "Very effective", "Extremely effective")



q33table <- svyby(~Q33A, ~raceparty, norcsvy, na.rm = T, svymean)
q33table <- full_join(q33table, svyby(~Q33B, ~raceparty, norcsvy, na.rm = T, svymean))
q33table <- full_join(q33table, svyby(~Q33C, ~raceparty, norcsvy, na.rm = T, svymean))
q33table <- full_join(q33table, svyby(~Q33D, ~raceparty, norcsvy, na.rm = T, svymean))

q33table <- q33table %>% melt() %>% separate(variable, sep = "\\([0-9]\\)\\s", into = c("Question", "Response")) 
q33spread <- spread(q33table, key = Response, value = value)

q33means <- subset(q33table, grepl("^Q", q33table$Question))


q33means$Response <- factor(q33means$Response, levels = c("Extremely effective", "Very effective", "Not at all effective", "Not too effective", "Somewhat effective"))
#indx <- q33means$Response %in% c("Somewhat effective", "Not too effective", "Not at all effective")
#q33means$value[indx] <- -1*q33means$value[indx]
indx <- q33means$Response == "Somewhat effective"
q33means$value[indx] <- 1/2*q33means$value[indx]

ggplot(data = subset(q33means, Question == "Q33A" & q33means$Response %in% highs), aes(x = raceparty, y = value, fill = Response)) + 
  geom_bar(stat="identity", position = "stack", width = 0.8, color = "black") + 
  geom_bar(data = subset(q33means, Question == "Q33A" & q33means$Response %in% lows), stat="identity", position = "stack", width = 0.8, aes(y = -value), color = "black") +
  theme_tufte() + coord_flip() + scale_fill_manual(breaks = efflevels, values = pal) + geom_hline(yintercept = 0, color = "white") +
  labs(title="Requiring all police officers to receive racial bias training", y="",x="") +
  theme(plot.title = element_text(size=14, hjust=0.5)) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(legend.position = "bottom") + scale_y_continuous(breaks=seq(-1, 1, .25))
ggsave(file = "output/q33a.pdf", width = 8, height = 2)

ggplot(data = subset(q33means, Question == "Q33B" & q33means$Response %in% highs), aes(x = raceparty, y = value, fill = Response)) + 
  geom_bar(stat="identity", position = "stack", width = 0.8, color = "black") + 
  geom_bar(data = subset(q33means, Question == "Q33B" & q33means$Response %in% lows), stat="identity", position = "stack", width = 0.8, aes(y = -value), color = "black") +
  theme_tufte() + coord_flip() + scale_fill_manual(breaks = efflevels, values = pal) + geom_hline(yintercept = 0, color = "white") +
  labs(title="Developing community policing programs in neighborhoods to help build relationships\nbetween the public and local police officers", y="",x="") +
  theme(plot.title = element_text(size=14, hjust=0.5)) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(legend.position = "bottom") + scale_y_continuous(breaks=seq(-1, 1, .25))
ggsave(file = "output/q33b.pdf", width = 8, height = 2)


ggplot(data = subset(q33means, Question == "Q33C" & q33means$Response %in% highs), aes(x = raceparty, y = value, fill = Response)) + 
  geom_bar(stat="identity", position = "stack", width = 0.8, color = "black") + 
  geom_bar(data = subset(q33means, Question == "Q33C" & q33means$Response %in% lows), stat="identity", position = "stack", width = 0.8, aes(y = -value), color = "black") +
  theme_tufte() + coord_flip() + scale_fill_manual(breaks = efflevels, values = pal) + geom_hline(yintercept = 0, color = "white") +
  labs(title="Offering incentives to members of the police force \nto live in the community where they work", y="",x="") +
  theme(plot.title = element_text(size=14, hjust=0.5)) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(legend.position = "bottom") + scale_y_continuous(breaks=seq(-1, 1, .25))
ggsave(file = "output/q33c.pdf", width = 8, height = 2)


ggplot(data = subset(q33means, Question == "Q33D" & q33means$Response %in% highs), aes(x = raceparty, y = value, fill = Response)) + 
  geom_bar(stat="identity", position = "stack", width = 0.8, color = "black") + 
  geom_bar(data = subset(q33means, Question == "Q33D" & q33means$Response %in% lows), stat="identity", position = "stack", width = 0.8, aes(y = -value), color = "black") +
  theme_tufte() + coord_flip() + scale_fill_manual(breaks = efflevels, values = pal) + geom_hline(yintercept = 0, color = "white") +
  labs(title="Requiring police departments to recruit additional qualified minority officers", y="",x="") +
  theme(plot.title = element_text(size=14, hjust=0.5)) +
  theme(axis.text.y = element_text(hjust=0)) +
  theme(legend.position = "bottom") + scale_y_continuous(breaks=seq(-1, 1, .25))
ggsave(file = "output/q33d.pdf", width = 8, height = 2)

